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1 . 


Introduction 


In this Report, research work done at Cambridge 
Hydrodynamics, Inc. under Contract No. NAS 1-15 3 72 with NASA 
Langley Research Center is summarized. More detailed 
expositions of the work described here are given in the 
publications cited in the References. 

In Sec . 2 , work done to extend the SALLY computer 
code (written by CHI with NASA contract support) to treat 
the propagation of wave packets is explained. 

In Sec. 3, work done to extend the SALLY computer 
code to treat nonlinear, nonparallel compressible flow over 
LFC wings is described. 



2 . Wave Packet Analysis for the SALLY Stability Analysis Code 


The SALLY computer code was developed for NASA 

Langley Research Center by CHI to be a 'black-box' stability 

analyzer for three-dimensional boundary layer flows over LFC 
1-3 

wings. Among the user options offered by SALLY are 

maximum growth rate determinations at (i) fixed frequency, 

(ii) fixed wavelength, (iii) fixed frequency, fixed wave- 
length, (iv) fixed wavelength, fixed propagation angle, 

(v) fixed frequency, fixed propagation angle. In the present 
work, SALLY has been extended to allow a wave packet analysis 
on an LFC wing. 

For steady boundary layer flow on an LFC wing, 
the equations of a wave packet are 


dx 

dt 


9 & 


( 1 ) 


dk 9u) 



( 2 ) 


where k is the wave vector of the packet and x is the 
current location on the wing. The required derivatives on 
the right side of ( 1 ) — ( 2 ) are computed using the adjoint 
eigenfunction of the Orr - Sommer f eld equation on the wing. 
Thus, if the Orr- Sommer f eld equation is written 


L (£,x, to (k,x) ) i/j = 0, (3) 

where ip is the eigenfunction, then differentiation with 
respect to k gives 
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(4) 


^ + L = 

31c 3a) 9k 9ic 

Taking the inner product of (4) with the adjoint x of the 
eigenfunction ip gives 


sS (x, §| *> 

Similarly, differentiation of (3) with respect 
to x gives 


so that 
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The modified SALLY code is being used by NASA LaRC 
personnel to study transition prediction using wave packets 
instead of maximum amplification rates. 
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3 . Compressible Flow Stability Analysis 


The SALLY stability analysis code has been extended 
to include the effects of compressibility. The equations for 
linearized disturbances to general three-dimensional com- 
pressible boundary layers are given in Refs. 2 and 4. Using 
the notation of Ref. 4, the equations may be written: 


where 



Y + A — ¥ + B¥ 
9y 


0 



( 8 ) 


(9) 


and A and B are 4x4 complex matrices whose elements 
are given below. Here a is the x-wave number, $ is the 
z-wavenumber , u and w are the amplitudes of x- and 

A 

z- velocity fluctuations, respectively, T is the amplitude 
of temperature fluctuations, and v is the amplitude of 
the normal velocity fluctuations. The elements of the A 
and B matrices are: 


A 11 

_ 1 dy , 
y dT 

(10) 

A 12 

a 2 / a i 

(11) 

A 13 

= ^ It (au'+ew’) 

(12) 

■ — i 
< 

= 0 

(13) 
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A 21 " a l a 3 a 4 


A 22 = “ a 5 + a 6 “ a 7 /a l 


A 23 = a 8 a 9 a l 


A 24 " 0 


A 31 a 10 


L 32 


a ll /a l 


A 33 a 12 


A 34 a 13 


A 41 = 0 


A 42 = 0 


A 43 a 14 


A 44 “ a 15 

B 11 = a l6 + a 4 a 2 /a l 
B 12 = a 17 + a 5 a 2 //a l 


Cl 4) 

(15) 

(16) 

(17) 

(18) 

(19) . 

( 20 ) 
( 21 ) 
( 22 ) 

(23) 

(24) 

(25) 

(26) 
(27) 
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( 28 ) 


B 13 a X8 - a 8 a 2 //a l 
B 14 = ° 

B 21 = a 7 a 4 //a l " a 19 ~ a l a 20 " a 4 a 6 

B 22 = a 7 a 5 //a l “ a 21 ” a l a 22 “ a 5 a 6 

B 23 = ~ a 7 a 8 //a l “ a 23 " a l a 24 + a 8 a 6 

B 24 = 0 

B 31 = “ a ll a 4 /a l 

B 32 = a 25 ” a ll a 5 //a l 

B 33 = a 26 + a ll a 8 /a l 

B 34 = 0 
B 41 = ° 

B 42 = a 27 

B 43 = a 28 

B 44 = a 29 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 
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(42). 


2 

a 1 = - iyM (aU+BW-oi). 

2 2 

a 2 = — (ot - j (l+2d) (a 2 +B 2 ) yM 2 (ctU+BW-oj) (43) 


a 3 

i. 

L 

(44) 

a 4 = 

- i 

(45) 

a 5 = 

T'/T 

(46) 

a 6 = 

% (2+d)YM^ (aU+BW-oi)— 3E t . 

■j e Jj y aT 

+ aU’+BW' + ^-(aU+ew-uj) J 

(47) 

a 7 = 

- iyM 2 (aU’+BW ) 

(4 8) 

a 8 

- ~(aU+BW-w) 

(49) 

II 

cn 

fd 

^ (2+d) ^ (aU+BW-w) 

(50) 

a io ~ 

2a(y-l)M 2 (aU’+BW )/(a 2 +B 2 ) 

(51) 

a ll = 

Rct o 

i — (y-l)M^ (aU+BW-w) 

(52) 

a 12 = 

2 dK 
K 3 t 

(53) 

a 13 = 

2a (y-l)M 2 (aW'-BU 1 )/(a 2 +B 2 ) 

(54) 
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(.55). 
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1 dy 
y dT 
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y dT 

iR 

yT 
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(aU+BW-co) - a 2 - B 2 


vT (aU'+BW)+i(a 2 +B 2 ) §£ T 


1 ? ? T 1 

+ ± i ( l+ 2 d) (cT+B ) 
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- | (l+2d) (a 2 +B 2 ) ^(aU+BW-oa) 


+ - ^ (aU"+BW") + - T' (ctU'+BW 1 ) 


y dT 


dT‘ 


a 19 = 0 


a 2o - - r [ ! It t ' + l (2+d) ^ 


T" (T ' ) 2 
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_ 1 r 2 Q 2 2, 0 , , . T ^ 1 dy 

a 22 L 01 ^ 3 ( 2 d ) o y dT 


+ | ( 2+d) ^ ~ (aU+BW-w) ] 


a 23 = I (aU'+BW) - (aU+BW-w ) 


(56) 
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(. 60 ) 
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(62) 

(63) 
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(69) 
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Here 

L = ^ + i|-(2+d) yM 2 (aU+BW-w) (71) 


and U , W are the boundary layer velocity components, T 
is the temperature, y is the viscosity coefficient, k is 
the thermal conductivity coefficient, a is the Prandtl 
number, y is the ratio of the specific heats, R is the 
Reynolds number, and d is the ratio of the second to the 
first viscosity coefficient. 

Equations (8)- (71) are to be solved subject to the 
boundary conditions of vanishing fluctuations at the 
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boundaries y = 0 , °°. 

The numerical method used to solve (_8).-(711 is 
explained in the Appendix. A new iteration procedure to 
solve the Chebyshev-spectral approximation to (8)— (71) has 
been developed that allows solution with up to 513 Chebyshev 
modes on the CDC Cyber computers at NASA LaRC. Some variants 
on the methods described in the Appendix are implemented in 
the compressible stability code. They are: 

(i) A variable resolution mapping in the 

y-direction is implemented. Early numerical 
experiments with the SALLY code indicated 
that the critical layer of compressible 
boundary layer flows lies much farther from 
the wall (wing) than in incompressible flows. 
In order to resolve adequately the flow region 
near this critical layer, a variable 
resolution mapping that gives extra weight 
to the critical layer is used. 

(ii) In order to minimize the computational expense 
due to resolution of oscillations of the 
boundary layer at large distances from the 
wing, where the boundary layer codes give 
constant unperturbed flow properties, the 
exact modes of oscillation of the boundary 
layer at large distances are determined by an 
algebraic eigenvalue analysis. The boundary 
conditions at 00 are then replaced by 
conditions at finite y that preclude the 
presence of growing oscillations at 00 . 
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Technical details of this and other matters 
concerning the compressible SALLY code will be given in a 
formal paper for publication in the open literature. This 
paper will be prepared upon completion of numerical 
calculations of the stability properties of various LFC 
wings using these compressible stability analysis codes. 
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4 . Nonlinear and Nonparallel Stability Analysis of 
Compressible Flows 

The equations for the nonlinear and nonparallel 
evolution of three-dimensional disturbances of three- 
dimensional compressible boundary layers were given in 
Ref. 2. Even with some notational simplifications, these 
equations require nearly 40 pages of text to write down. 

Their direct coding on a computer would be extremely 
laborious and likely to be subject to serious error. 
Fortunately, a much easier way to solve the equations re- 
quired Cor, rather, to develop a theory that is nearly the 
same as that of Ref. 2 that leads to much simpler equations 
than in Ref. 2) has been developed. 

The key idea is to apply the transform methods 

5 

developed by CHI for solution of the Navier-Stokes equations . 
In this Section, a simple example of this idea is given 
postponing full discussion of the technical details of its 
application to compressible flow until a formal paper on the 
results is prepared for publication. 

Consider the inviscid Burgers' equation 


3u(x,t) , /,, . \ 3u(x,t) _ q 

3x ~ U 


(72) 


with periodic boundary conditions u(x+27r,t) = u(x,t)'. Let 
us consider perturbations to the steady solution u(x,t) = U 
of the form 


u = U + e Re (u^ ( t) e^ kx ) 

+ e 2 Re(u 2 a)e 2ikx +u 2 e ikx ) 

*5 ”5 *i Iry 

+ s Refuse x +»**) + ••• 


(73) 
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When nonlinear stability theory is applied to this prohlem, 
it is necessary to compute the following terms: 

iUu l 

2iUu^ + i ( ) 2 (74). 

3 2 1 * 

iUu^ + iu 2 (u^) , etc. 

For the present example, the computation of these terms is 
obviously straightforward and presents no unusual 
difficulties. However, the equations of compressible 
stability analysis are considerably more complicated and the 
calculation of the required terms presents no trivial 
challenge to the computer. 

The new method involves the evaluation of the non- 
linear term u 3u/3x using collocation methods at discrete 
points x^ and using the explicit form (73), neglecting 

all terms not shown explicitly. Inverse fast Fourier trans- 
formation then allows the computation of wavenumbers 1 , 2 
of the nonlinear term (still as a function of £^). If the 
results of this pseudospectral computation are denoted N(£j), 
then the required nonlinear interaction terms in (74) are 

3 

computed as follows. If terms through order e are required, 
then N should be evaluated using 4 arbitrary values of 
e j . The various coefficients are determined by solving for 
Nq , , N 2 , N 3 from 

N(e.) = N + N. £ . + N £? + N £7 + ••• . (75) 

3 u ±3 ^3 J 3 
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Choosing = z/2? for some small e > 0 allows the 
efficient computation of the nonlinear interaction terms 
V N^, N 2 , using the formulae for Richardson extra- 
polation. ^ 


Using these ideas, it has been possible to 
code the nonlinear and nonparallel flow terms for com- 
pressible boundary layers very efficiently, with minimal 
chance for programming error. The computer codes involve 
just the programming of the nonlinear terms of the com- 
pressible Navier-Stokes equations by pseudospectral 
Ccollocation) methods and then a final Richardson extra- 
polation to extract the required information. These new 
subroutines are being installed in the standard SALLY code 
during FY8Q . 

Similar ideas to those outlined above should 
revolutionize numerical methods for nonlinear, nonparallel 
stability theory in more general problems. 
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Appendix 


1. INTRODUCTION 


In this appendix some new techniques . are introduced that permit 
the efficient application of spectral methods to solve proo- 
lems in (nearly) arbitrary geometries. The resulting methods 
are a viable alternative to finite difference and finite 
element methods for these problems. Spectral methods should 
be particularly attractive for problems in several space 
dimensions in which high accuracy is required. 

Spectral methods are based on representing the solution 
to a problem as a truncated series of smooth functions of the 
independent variables. Whereas finite element methods are 
based on expansions in local basis functions, spectral methods 
are based on expansions in global functions. Spectral methods 
are the extension of the standard technique of separation of 
variables to the solution of arbitrarily complicated problems. 


Spectral methods are best introduced for the 
simple one-dimensional heat equation. Consider the mixed 
initial-boundary value problem 


3u(x,t)_ Tr 3 2 u(x,t) 
' 3x 2 


(0<X<TT , t>0) 


( 1 . 1 ) 


U(0,t) = U(7T,t) = 0 (t>_0) 

u(x,0) = f (x) ((Kx<rr) 

The solution to this problem is 

OO 

u(x,t) - Fa (t) sin nx 
n 

n=i 

a (t) = fe~ Kn 
n n 

where 

2 r* 

f - — f f(x)sinnx dx 
n it . J 

0 


( 1 . 2 ) 

(1.3) 


(1.4) 


(1.5) 


( 1 . 6 ) 
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are the coefficients of the Fourier sign series expansion of 
f (x) . 

A spectral approximation to (1.1) - (1.3) is gotten by 
simply truncating (1.4) to 

N 

u M (x,t) = la (t) sin nx (1.7) 

N n=l n 


and replacing (1.5) by the evolution equation 


da 

n 

dt 


Kn 2 a 


n 


(n=l , 


,N) 


( 1 . 8 ) 


with the initial conditions a (0) = f (n=l,...,N). 

n n 

The spectral approximation (1.7) - (1.8) to (1.1) ~ 

(1.3) is an exceedingly good approximation for any time t 
greater than zero as N+ 03 . In fact, the error u(x,t) - 
u^(x,t) satisfies 

oc 2 2 

u(x,t) - u (x,t) = £ f e Kn t sin nx = 0(e KN 

N n=N+l n 

(N+“) (1.9) 

for any t>0. In contrast to (1.9), finite difference 
approximations to the heat equation using N grid points in 
x lead to errors that decay only algebraically with N as 
N+°°. Furthermore, this spectral method for the solution of 
the heat equation is efficiently implementable by the fast 
Fourier transform (FFT) in 0(N £og N) operations. 

There are several significant difficulties in extending 
the simple spectral method employed for (1.1) - (1.3) to 

more general problems. These difficulties and their solutions 
will be discussed in the following sections. Some further 
details are given in the author's monograph [1) . In Sec. 2, 
the difficulty caused by nontrivial boundary conditions 
is discussed. In Sec. 3, the difficulty of treating nonlinear 
and nonconstant coefficient terms is discussed. Then, in 
Sec. 4, the properties of spectral methods for problems in 
simple geometries are summarized. In Sec. 5, the extension of 



spectral methods to problems in complicated geometries is 
explained. In Sec. 6, a new technique for the efficient 

solution of spectral equations that arise in complicated 
geometries is given. Some representative test problems are 
discussed in Sec. 7. Then, in Sec. 8, a summary of results 
is given together with a glimpse of some other new develop- 
ments in spectral methods that should find wide application. 


2. ' BOUNDARY CONDITIONS 

The Fourier series (1.4) converges fast if u(x,t) is 
infinitely differentiable and u(x,t) satisfies the boundary 
conditions 


3 2n u(x,t) 


3x 


2n 


0 (x = 0 ,it) 


( 2 . 1 ) 


for all nonnegative integers n. Under these conditions, the 
error after N terms 

N 

e N (x,t) = u(x,t) - l a n (t) sin nx 


goes to zero uniformly in x faster than any power of 1/N 
as N -* °o . On the other hand, if u(x,t) is not infinitely 
differentiable or if any of the conditions (2.1) is 
violated, then s N (x,t) = 0 (1/N P ) as N -* <=° for some 
finite p. For example, 

1 = l (-1) — y— i — ( 0<x<ir ) , (2.2) 

n=0 zn+x 

but the error incurred by truncating after N terms is of 
order 1/N for any fixed x,0<x<tt. Furthermore, the 
convergence of (2.2) is not uniform in x; (2.2) exhibits 
Gibbs’ phenomenon, namely 

s n (C/N) = 0(1) (N+°°, £ fixed). 

For any fixed N, there are points x at which the error 
after N terms of (2.2) is not small. The poor convergence 
of (2.2) is due to the violation of (2.1) for n = 0. 



More generally, most eigenfunction expansions of a 
function f (x) converge faster than algebraically (i.e. 
the error incurred by truncating after N terms goes to zero 
faster than any finite power of 1/N as N-»-<») cnly if f (x) 
is infinitely differentiable and f (x) satisfies an infinite 
number of special boundary conditions. For example, the 
Fourier-Bessel expansion 

00 

f(x) = l a n J Q U n x) (0<x<l) 


where A is the nth smallest root of J_(X) = 0, converges 
n u 

faster than algebraically only if f is infinitely different- 
iable and 




1 

x 


d 

dx 


x 



f (x) 


0 at x = 1 


(2.3) 


for k = 0,1,2,... . 

When a spectral expansion converges only algebraically 
fast, spectral methods based on these eigenfunction expansions 
can not offer significant advantages over more conventional 
(finite difference, finite element) methods. Eigenfunction 
expansions of this kind should not normally be used unless 
the boundary conditions of the problem imply all the extra 
boundary constraints like (2.1) or (2.3). For example, 
if periodic boundary conditions are compatible with the 
differential equation to be solved, complex Fourier series 
are suitable to develop efficient spectral approximations. 

In the development of spectral methods for general 
problems, it is important that the rate of convergence of the 
eigenfunction expansion being used not depend on special 
properties of the eigenfunctions, like boundary conditions, 
but rather depend only on the smoothness of the function 
being expanded. Of course, if the solution to the problem 
being solved is not smooth, one should not expect errors 
that decrease faster than algebraically with 1/N when 
global eigenfunction expansions are used. Faster than 
algebraic rates of convergence may be achieved for these 
problems by either patching the solution at discontinuities 
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(see Sec. 5) or pre-and post-processing of the solution 
(see [2] ) . 

There is an easy way to ensure that the rate of converg- 
ence of a spectral expansion of a function f (x) depends only 
on the smoothness of f (x) , not its boundary properties. The 
idea is to expand in terms of suitable classes of orthogonal 
polynomials, including Chebyshev and Legendre polynomials for 
all those problems in which constraints like (2.1) and 
(2.3) are unrealistic. These polynomial expansions avoid all 
difficulties associated with the Gibbs phenomenon provided 
the solution f (x) is smooth. 

From the mathematical point of view, the classical 
orthogonal polynomials are eigenfunctions of singular Sturm- 
Liouville problems. It is not hard to show (see [1] for the 
details) that expansions using eigenfunctions of such sing- 
ular Sturm-Liouville problems converge at a rate that depends 
only on the smoothness of f (x) , in contrast to eigenfunction 
expansions based on nonsingular Sturm-Liouville problems that 
lead to additional boundary constraints like (2.1) on f(x). 

These results for orthogonal polynomial expansions are 
easily demonstrated in the case of Chebyshev polynomial 
expansions. The nth degree Chebyshev polynomial T n (x) is 
defined by 

T n (cos 0 ) = cos n 0 . (2.4) 

Therefore, if 

00 

f(x) = l a T (x) ->(2.5) 

n=0 n n 

Then 

00 

g(0) = f (cos 0 ) = £ a cosn0 (2.6) 

n=0 n 


Thus, the Chebyshev polynomial expansion coefficients a 

n 

of f (x) are just the Fourier cosine expansion coefficients 
of the even, periodic function g(6) . A simple integration 
by parts argument then shows that 



(n ■+ °°) 


a 0 
n 

provided g(0) (or, equivalently, f (x) ) has. p continuous 
derivatives. Since 

|f(x) - l a T (x) | < l la | (|x|.<l), 

n=0 n n n=N+l n 


it follows that the rate of convergence of (2.5) is faster 
than algebraic if f is smooth. 

In summary, spectral expansions should be made using 
series of orthogonal polynomials unless the boundary conditions 
of the problem are fully compatible with some other class of 
eigenfunctions. In practice, Chebyshev and Legendre polynomial 
expansions are recommended for most applications, supplemented 
by Fourier series and surface harmonic series when boundary 
conditions permit. 

3. NONLINEAR AND NONCONSTANT COEFFICIENT PROBLEMS 

Another difficulty with general kinds of spectral 
methods is their application to problems with nonlinear and 
nonconstant coefficient terms. Before explaining the solution 
to this problem, an illustration of the difficulty is helpful. 
Consider the partial differential equation 


ff = N(u,u) + L u ;(3.1) 

* y 

where u = u(x,t) and N is a bilinear (nonlinear) operator 
that involves only spatial derivatives and L is a linear 
operator that involves only spatial derivatives. The 
operators N and L may depend on both x and t. A 
spectral method for the solution of (3.1) is obtained by 
seeking the solution as a finite spectral expansion: 

N 

u(x,t) = l a (t) if# (x) (3.2) 

n=l n n 
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where it is assumed for now that ^ n (x) (l<n<«>) are a complete 

set of orthogonal functions. Defining the re-expansion coefficients 

c and d by 

nmp run J 


= £ c (t) ip 

m p nmp n 


(3.3) 


(-(<!') = I d (t) ip 


(3.4) 


m 


n=l 


nm 


n 


and equating coefficients of ij; (x) (n=l,...,N) in (3.1), 

gives 


da 

n 

dt 


N 

l 

m=l 


N 

I 

p=l 


c (t) a (t) a (t) 
nmp m p 


N 

+ I 

m=l 


nm 


(t) a (t) 
m 


(n=l , . . . ,N) (3.5) 


Eqs . (3.5) are the spectral evolution equations for 

the solution of (3.1). They have one very serious drawback. 

In general c nmp and d nm are nonzero for typical n,m,p 

so that evaluation of da /dt from (3.5) for all 

3 n 

n=l,...,N requires 0 (N ) arithmetic operations for the 

2 

bilinear term and 0 (N ) operations for the linear term. 

3 

Thus# solution of (3.5) requires order N operations per 
time step. Since operational spectral calculations now 
involve N £ 10^, the computational cost of the direct 
solution of (3.5) is prohibitive (even if only linear terms 
are present) . 

The problem here is one of computational complexity. 

Finite difference methods for the solution of (3.1) on N 

grid points may require only order N operations per time 

3 

step. If the spectral method really requires order N 
operations per time step it can not compete when N is 
large . 

Another example illustrating the computational complexity 
of spectral methods is given by the nonlinear diffusion 
equation . 


9u (x, t) 
8t 


u 3 u 


3x' 


(x,t) . 


(3.6) 
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I 


Expanding the solution as 

N 

u(x,t) = l a (t)i|j n (x) (3.7) 

n=l 


in terms of the orthonormal functions ^(x) ' gives 

N 

I 


for n = 1,...,N. These evolution equations for {a n (t)} 
have an exponential degree of computational complexity as 
they are expressed as an integral functional of {a^ (t) } . 
The solution to the problem of computational complexity 

is to use the author's transform methods. This technique may 
be illustrated for a pseudospectral (or collocation) approx- 
imation to (3.6) [3], First, N suitable collocation points 

x ,x , ...,x lying within the computational domain are intro- 
■L Z N 

duced. Then, the approximate solution (3.7) is forced to 

satisfy the partial differential equation (3.6) (or its 

boundarv conditions), exactlv at these discrete ooJ.nts at 

every time t. More specif ically , the following three steps 

are done at each time step t: 

(i) Determine N coefficients a (t) (n-l,...N) so 

n 

that 

N 

u (x . , t ) = l a. (t)ty (x.) (j=l, ,N) . (3.9) 

n=l 3 


a (t) ip (x) 
m m 


N 

l 

p=l 


a ' 
P P 


(x) dx 


da 

-gr = / V x) exp 


(ii) Evaluate u (x.,t) by 

xx j 


u (x . , t) 
xx 3 


N 

l a (t)^» " (x ) (j=l / ...,N). 

n=l n n J 


(3.10) 


(iii) Finally, evaluate 


3u (Xj , t) /3t 


3u (x_. , t) 


3t 


u (x . ,t) 

e J u (x . , t) ( j=l , . . 
xx j J 


by 

,N) 


(3.11) 
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and march forward to the next time step. 

The idea of the pseudospectral transform method can be 
restated as follows: Transform freely between physical (x j ) 

and spectral (a n ) representations, evaluating each term in 
whatever representation that term is most accurately, and 
simply, evaluated. Thus, in (3.11), e U is evaluated in 
the physical representation while u is compared in the 
spectral representation by (3.9) because it is most accurately done 
there. 

It should be apparent to the reader that pseudospectral 
transform methods can be applied to any problem that can be 
treated by finite difference methods regardless of the 
technical complexity of nonlinear and nonconstant coefficient 
terms . 

What is the computational complexity of pseudospectral 
transform methods? There are at least three 

aspects to this question: (i) the computational complexity 

of differentiation, integration, etc. of spectral series; 

(ii) the computational complexity of transforming back and 
forth between physical and spectral representations; and 

(iii) the computational complexity of solving the resulting 
equations for the spectral coef f icients . 

For the expressions of interest, computation of deriva- 
tives of an N term spectral expansion requires order N 
arithmetical operations. For the Fourier series (1.7), this 
fact is obvious: 

* N 

d v 

j— > a sm nx 

dx n 

n=l 


N 


n=l 


n a cos nx 
n 



N 


l 

n=l 


a sin nx = - 
n 


N 


l 

n=l 


2 

n a sin nx . 
n 


For the Chebyshev polynomial expansion (2.5), the comput- 
ational complexity of differentiation is a little less apparent . 

Since T (cos 0 ) - cos n 0 , 
n 


T A+i <*> 

n+1 


T 1 , (x) 
n-1 

n-1 


— T (x) (n>0) 

c n ~ 
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I 


where c n = 2, c =1 (n>l) and T ' = T' = 0. Therefore, 

On— u — 1 

if 


, N N 

-=— y a T (x) = 7 b T (x) , 

dx n n L n " 


n=0 


n n 


then 


N 

2 y a T ' (x) 

nil n n 


v k r t a+i T A-ii 

n | 0 C n b n L n+l n-lj 


N+l 

y [c -i b , -b _ ] T ' (x)/n 
n-1 n~l n+l n 

n=l 


Equating coefficients of T^(x) for n = 1,...,N+1 gives 
the recurrence relation 

c , b , - b . , = 2na (l<n<N) 

(3.12) 


given a^ requires only 
Similar recurrence relations 
can be obtained for differentiation of spectral series based 
on other sets of orthogonal polynomials and functions. 

The computational complexity of transforming between 
spectral and physical space has several interesting aspects. 
The problem is: how much computational work is necessary 

to evaluate 

N 

u. = y a * (x.) (j=l,...,N) (3.13) 

D n =l n n 3 


n-± n-r n+± n 

b = 0 (n>N) 

n — 

The solution of (3.12) for b^ 
order N arithmetic operations. 


given {a^} and, inversely, how much work is necessary to 

compute the expansion coefficients {a^} given {u^ } ? 

It is obvious that (3.13) can be evaluated for {u.} in 
2 3 3 

0 (N ) operations while it can be solved for a in 0 (N ) 
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operations. However, these estimates are much too pessimistic; 
for many important expansions the operation count to perform 
the transform (3.13) and its inverse is no larger than 
0 (N ( log N jP) with p £ 2 . 

In the case of Fourier series, the transform (3.13) and 
its inverse can be computed in 0 (N 5og2 N) operations of 
N = 2 P using the fast Fourier transform. However, most of 
the computational efficiency of transform methods comes not 
from the FFT but from the separability of multi-dimensional 
transforms. Thus, a three-dimensional discrete Fourier trans- 
form can be expressed as three one-dimensional Fourier trans- 
forms 


J-l 

l 

j=o 


K— 1 L-l 

l l a(j,k,£) exp[2TTi(i^ + ~r + 4jr) J 
k=0 1=0 


(3.14) 


e 2iTijm/J e 27rikn/K 

j-0 k=0 


... 0 . 2tt i£p/L 

l a (j ,k, £,) e 

5,-0 


2 

The left side of (3.14) requires roughly ( JKL ) 
operations to evaluate at all the points 0<m J, 0<n<K , 0<p<L . 
On the other hand, even without the FFT the right side of 
(3.14) requires only about (JKL) (J+K+L) operations to 
evaluate at all the points. When the FFT is applied to the 
one-dimensional transforms on the right side of (3.14), the 
number of operations necessary to evaluate (3.14) is reduced 
further to (JKL) (&og 2 J + ^ 0< ?2 K + ^ J,K,L are 

powers of 2. 

For example, the spectral code CENTICUBE solves the 
Navier-Stokes equations for three-dimensional incompressible 
flow with periodic boundary conditions; three-dimensional 
Fourier series with resolutions up to 128 * 128 * 128 are 
used to represent the flow. The equations for the Fourier 
components u(£,t) of the velocity field involve several 
convolution sums of the form 


l u ( p , t ) u (k-p, t) . 

|p| <64 
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It 'is not hard to show that direct evaluation of all the 

required convolution sums on the CRAY-1 computer would require 

5 

(using optimized code) about 5 xlO s of computer time at 

each simulation time step. On the other hand, the CENTICUBE 

code executes on the CRAY-1 in less than 20s per time step! 

This speedup by a factor 2.5 *10^ is broken down roughly 

as follows: a factor 2 for using the FFT instead of an 

optimized matrix multiply to perform a one-dimensional trans- 

4 

form and a factor 10 to perform the transforms as a seq- 
uence of one-dimensional transforms as in (3.14), 

More general transforms can also be implemented efficient- 
ly. The author has recently shown that 'fast' trans- 
forms of N-term series of Legendre polynomials, surface 
harmonics, Bessel functions, Jacobi polynomials, Gegenbauer 

polynomials, etc. cap be accomplished in 
2 

0(N(£og N) / £og Jlog N) arithmetic operations while trans- 
forms of series of Hermite and Laguerre polynomials can be 

2 

accomplished in 0(N(£og N) ) operations. At the present 
time, these transforms are of only academic interest - as 
discussed above the speed of transform methods for problems 
of realistic size owes to the formulation of the problem in 
terms of separable transforms not the existence of a fast 
transform. 

The third question regarding computational complexity of 
spectral transform methods concerns the complexity of solving 
the equations for spectral coefficients. In the case of 
initial value problems solved by explicit time-stepping 
methods, the answer is provided by the transform methods 
discussed above: 0(N(£og N)^*) operations are required per 

time step. The answer to the question is more complicated for 
the solution of boundary value problems or implicit time- 
stepping methods fo.r initial value problems. 

Spectral approximations to general boundary value problems 
lead to full N x N matrix equations for the N expansion 

coefficients a . It would seem that solution of these 
n 3 

equations requires 0(N ) arithmetic operations while 

2 

storage of the matrix requires 0 (N ) memory locations. 

6 

Since typical problems now involve N ~10 , the direct 
solution (or even the direct formulation) of such problems 
is clearly unworkable now. The solution to this problem is 
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given in Sec. 6. The solution to the spectral equations 
requires essentially only 0(N(£og N)^) (p^_2) operations 

and only 0(N) memory locations. Solution of spectral 
equations, even though they lead to formally full matrix 
problems, can be accomplished in little more work and memory 
than required to solve the simplest finite difference equations 1 

4. TIME DIFFERENCING AND BOUNDARY LAYERS 

Spectral methods based on the classical orthogonal poly- 
nomials have another feature that is very attractive for some 
kinds of problems but leads to difficulties with others. This 
feature is very high resolution near the boundaries. 

For example, the collocation points for Chebyshev polynomial 
pseudospectral methods for problems on -1<x<1 are usually chosen 

to be x. = cos "irj/N { j = 0 , 1 , . . . , N) . The collocation points 

3 .22 

x^ and are within about ir /2N of the boundary 

points Xq and x^ , respectively, so that the boundary 

resolution is Ax = 0{1/N^). This leads to extremely 

good resolution properties of spectral methods for boundary 

layer problems (see [1] and [4]) . While resolution of a 

problem with a boundary layer of thickness e<<l would 

require 0(l/e) uniformly spaced grid points, it requires 

only 0(l//e) terms in the Chebyshev spectral series. [Non- 

uniform grids should, of course, be used in many of these 

problems. They can also be implemented efficiently in 

spectral methods using coordinate transformations.] 

The high boundary resolution of spectral methods is not 

directly useful when problems without boundary layers are to 

be solved. For example, consider the one-dimensional wave 

equation 


u t + u x 

= 0 

*(-l<x<l , t>0 ) 

(4.1) 

u(-l,t) 

= f (t) 

(t>0) 

(4.2) 

u(x,0) = 

: g(x) 

(-1<x<1) . 

(4.3) 


Using a Chebyshev polynomial spectral representation 
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N 

u(x,t) = l a (t)T (x) (4.4) 

n=0 


where T (x) = 
n 

the solution of 


da 
n 

dt 



cos(n cos 
(4.1) are 

N 

l pa 

p=0 

p+n odd 


x) , the spectral-tau equations for 

[ 1 ] : 

(0<n<N-l) (4.5) 

p 


N 

l (-1) n a (t) = f(t) 
n=0 n 


(4.6) 


where c n = 2, c ~ 1 (n>l) . The numerical solution of 
(4.5) - (4.6) in time can be done using any absolutely stable 
time integration method for ordinary differential equations, 
such as one of the Adams methods [5] . For an explicit Adams 
method (Adams-Bashforth method) , absolute stability in 
solution of (4.5) - (4.6) requires that 


At = 0 (-—) . 
IT 


(4.7) 


This result may be verified from (4.5) - (4.6) .using 

Gerschgorin ' s theorem on the distribution of eigenvalues 

or, more physically, from the classical explicit stability 

2 

condition At = 0(Ax) with Ax = 0(1/N ). Specifically, 
the first-order Adams-Bashforth method (Euler's method) is 
stable for the solution of (4.5) - (4.6) provided that 

a t< . (4.8) 

tr 


The time step restrictions (4.7) - (4.8) are quali- 
tatively more severe than those encountered by standard 
finite difference methods for (4.1) - (4.3). The solution 
to (4.1) - (4.3) does not exhibit boundary layer structure 
(unless g(x) or f(t) have non-uniform variation) so a 
uniform grid may be employed giving the stability criterion 
At = 0(1/N) using N grid points. The high boundary 
resolution of the spectral scheme that leads to the more 
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stringent requirements (4.7) or (4.8) may seem wasted in 
this problem. In fact, this high boundary resolution is not 
completely useless; while high-order accurate stable finite 
difference schemes for solution of (4.1) - (4.3) on a 
uniform grid are complicated and require a number of 
spurious numerical boundary conditions (see, for example, [6]), 
the infinite-order accurate spectral scheme (4.5) - (4.6) 
is quite straightforward and requires no spurious boundary 
conditions to be applied. However, it is also nice to know 
that the stiffness of the spectral equations can be easily 
circumvented and time step restrictions like those of finite 
difference schemes can be easily obtained. 

At the present time, there are three alternative ways 
to avoid the severe time step restrictions of spectral methods 
in problems that do not exhibit strong boundary-layer structure. 

First, a semi-implicit scheme [1] may be employed to relax 
the time step restrictions (4.7) - (4.8) to At = 0(1/N) 
as for finite difference schemes. The idea here is to 
treat implicitly just those parts of the problem, to wit, the 
boundary regions, that lead to the stiffness of the spectral 
equations . 

Second, it is possible to formulate explicit uncondit- 
ionally stable time differencing schemes for spectral methods (Turkei & 
Gottlieb, to be published) . Domain of dependence arguments may be used 
to demonstrate the existence of conditional stability bounds for 
explicit finite difference methods. These bounds can be 
avoided in explicit spectral methods because global data is 
involved at each time step to march the solution forward in 
time. A simple example of an unconditionally stable explicit 
spectral method may be given for 

= -iku(k,t) (4.9) 


which is the equation for a Fourier mode of the solution 
to 3u/9t + 3u/3x = 0 with periodic boundary conditions. 
Leapfrog time differencing of (4.9) gives 


u (k, t+At) -u (k, t-At) 
2At 


iku (k , t ) 


(4.10) 
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which is second-order accurate in 


stability bound 


kAt I <1. 


At and has the explicit 
However the modified scheme 


u (k , t+At) -u (k , t- At ) _ 
2At 


i 5iSkAt uCk(t) 


<4.111 


is still second-order accurate but it is unconditionally 
stable since |sinkAt|_<l for all kAt. [By an accidental 
cancellation (4.11) happens to give the exact solution 
e -ikAt Qf (4>9) for all kAt.] 

The third method to stabilize time integrations of 
spectral methods is to use a full implicit time integration 
method. With N degrees of freedom used to resolve the 
spatial dimensions, full implicit schemes give full N x N 
sets of linear equations to solve at each time step. An 
efficient method to solve these equations is presented in 
Section 6. 

5. FORMULATION OF SPECTRAL METHODS IN COMPLEX GEOMETRIES 

For most problems in complex geometries, it is both 
inefficient and disadvantageous to seek special sets of 
spectral expansion functions tuned to the details of the 
geometry. First, determination of such a set of special 
functions is itself a computationally difficult problem in 
the complex geometry that must be repeated for every new 
qeometry. Second, unless very special care is taken, the 
resulting soectral series will not be guaranteed of fast 
convergence properties for the problem of interest. Third, 
much storage will be required for the expansion functions them- 
selves and fast, separable, transforms between physical space 
and transform space will not normally be available. 

There are two very general and powerful methods for the 
formulation of spectral methods in complex geometries that 
appear to preserve all the nice properties of spectral methods 
in simpler geometries, namely, mapping and patching. 

The idea of mapping is to transform the complex geometry 
into a simpler one by a smooth transformation. For example, 
the annular region 
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(5.1) 


( 0) < r < f 2 (0) 


0 < 0 < 2tt 


where (r,0) are polar coordinates is tranformed into the 
rectangle 


-1 < z < 1 


0 < 0 < 2tt 


(5.2) 


by the simple stretching transformation 
r-f i ( 0 ) 

2=2 f 2 (e)-f 1 (e) " 1 * 


(5.3) 


In the mapped coordinates (5.2), a spectral expansion using 
Chebyshev polynomials in z and Fourier series in 0 is 
appropriate (because the solution to the problem must be 
periodic in 0) . The complication of using a coordinate 
transformation appears in the coefficients of the differential 
equation in the transformed domain. For example, derivatives 
are transformed according to 


9u , 


2 


3u 



9r 

0 

f 2 ( 9 ) _f i 

(6) 

9 z 

0 


9u | 


9u | 

(z(f 2 

i 

- f ] 


9u 

90 | 

r 

90 

1 z 


f 2~ 

- £ 1 

3z J 


(5.4) 


(5.5) 


The complication of the equation due to transformation 
causes no essential difficulty in the implementation of 
explicit time stepping schemes for initial-value problems 
because transform methods are still applicable. For boundary 
value problems and full implicit treatment of initial value 
problems, it is essential that the full matrix equations be 
solved by the fast iteration procedures introduced in the 
next Section. 
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The simple stretching transformation (5.3) can only 
be applied if the boundaries r = f^ and r = f 2 are single- 
valued functions of 0. More general boundaries require more 
sophisticated mappings. In two dimensions, conformal mapping 
is sometimes useful. Also ,pseudo-Lagrangian marker particles 
may be used to define a coordinate frame. Many other sources 
of suitable transformations are possible (see, for example, 

[7]) . 

The idea of patching is to subdivide a complicated 
region into a number of simpler regions, make spectral expan- 
sions in each of the simpler regions, and then solve the 
problem in the complicated region by applying suitable 
continuity conditions across the artificial internal bound- 
aries. For example, consider the Poisson equation 

V 2 u = f 


in the L-shaped domain 

0 £ x £ 1 , 0 < y < 1 
-1 < x £ 0,0 1 y 1 2. 

The domain is subdivided into three domains: 


(5.6) 


1 • 0 < x < 1, 0 < y £ 1 (5.7a) 

II: -1 £ x £ 0,0 £ y £ 1 (5.7b) 

III: -1 £ x £ 0,1 £ y £ 2 (5.7c) 


In each of these regions, u(x,y) is represented as a double 
Chebyshev series : 


u I (x,y) = l 

u n ( x ,y) = I 


$ T m (2x-l) Tn (2y-l) 

II 

l a T (2x+l)T ( 2y-l ) 
L mn m n J 


(5.8a) 

(5.8b) 
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(5.8c) 


III 


u III (x ' y) = ? ^ a mn 


T (2x+l)T (2y-3) 
m n 


Across the internal artificial boundaries, continuity of u 
and the flux of u is imposed. For example, at the interface 
between regions I and II, it is necessary that 

u I (0,y) = u XI (0,y) (0 £ y £ 1) (5.9a) 



(0,y) 


9u 


II 


9 x 


(0,y) 


(0 < y < 1) . 


(5.9b) 


The system of equations in regions I, II, III together with 
the continuity conditions of the form (5.9) gives a 
spectrally patched solution to the Poisson equation in the 
L-shaped domain (5.6). 

The advantage of spectral methods over more convent- 
ional methods for patched problems is that the spectral 
schemes require only physical continuity conditions at the 
internal interfaces. On the other hand, finite difference 
methods require spurious boundary conditions at interfaces 
whenever the order of the numerical scheme is higher than 
that of the differential equation to be solved. 

6. SPECTRAL ITERATION METHOD 


Consider the solution of a general linear differential 
equation Lu = f. (Extensions to nonlinear problems will be 
discussed later.) Let an N-term spectral approximation to 
this problem be given by 


L sp ^ f N 


( 6 . 1 ) 


where f„ T is a suitable N-term approximation to f. As 
N 

mentioned several times earlier the matrix representation of 

(6.1) is generally a full NxN matrix so that direct 

solution of (6.1) by Gauss elimination methods would require 
2 

order N storage (for the matrix representation of L ) 

3 ^p 

and order N arithmetic operations. 
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In this Section, a method will he described that permits 
solution of (6.1) using order N storage locations with 
the number of arithmetic operations of order the larger of 
N £og N and the number of operations required to solve 
Lu = f by a first-order finite-difference method. The 
important conclusion is that spectral methods for general 
problems in general geometries can be implemented efficiently 
with operation costs and storage not much larger than that of 
the simplest finite-difference approximation to the problem 
with the same number of degrees of freedom . Since spectral 
methods require much less degrees of freedom to achieve a 
given accuracy or achieve much higher accuracy for a given 
number of degrees of freedom than required by finite-order 
finite-difference approximations, important computational 
efficiencies result 'from the new method. 

The idea of the iteration method is as follows: Suppose 

that it is possible to construct an approximation L to the 

ap 

spectral operator L that has the following properties: 


(i) L ^ as a s P arse matrix representation so that 

0 (N) 


it can be represented using only 

(ii) 

equation 


storage locations; 

L is efficiently invertible in the sense that the 
ap 


L U vr = f x-r 

ap N N 


( 6 . 2 ) 


is soluble as efficiently as a first-order finite-difference 
approximation to the problem. 


(iii) 


L approximates 
ap 


0 < m < II L*" 1 L 

— * 1 ap sp 


sp 
< M < 00 


in the sense that 


(6.3) 


for suitable constants m, M as N -»■ <*. Roughly speaking, 

(6.3) requires that the eigenvalues of L ^ L be bounded 

ap sp 

from above and below as N -*■ °°. Methods to construct suitable 


operator approximations L 


ap 


will be discussed below, 
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that 


There are several iteration procedures using L 


ap 


converge efficiently to the solution of (6.1). Three of these 

procedures are described below: 

Richardson (Jacobi) iteration 

Consider the iteration scheme 


L u* n+1) = L u* n) - a(L u^ n) - f ) 
ap N ap N sp N N 


(6.4) 


If 


0 < a < g 


(6.5) 


then 


u (n) ->- u (n-^-oo) 
N N ’ 


( 6 . 6 ) 


where u„ satisfies (6.1) : L u T = ft. The proof is 

N (n) (n) sp N ^ 

elementary: If e ' = “ U N i- s error / then 

e (n+l) _ (n) _ aL -l L e (n) 

ap sp 


Therefore, noting (6.3), 


| | e (n+1) | | £ max ( 1 1-am | , | l-aM| ) | | e | | 


(6.7) 


so 


(n) 


||e wl/ ||^0 as n 00 if (6.5) holds 
The optimal rate of convergence of Richardson’s iteration 
is normally achieved by choosing a to minimize the factor 
max ( | 1-am | , | 1-aM | ) appearing in (6.7). This gives 


opt 


M+m 


( 6 . 8 ) 


so 


(n+l) 

WT7 


M-m 
— M+m' 


(6.9) 
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Since, as shown below, it is possible to find L such that M <_ 2 . 5 

and m > 1 for nearly arbitrary spectral operators L , it 

sp 

follows that Richardson's method decreases the norm of the 
error in the solution of (6.1) by at least a factor 
(5/2+1) / (5/2-1 ) = 7/3 at each iteration. (Here ci 0 p-t = 4/7.) 

The accuracy of a typical initial approximation to u N is 
improved by a factor 10^ after about 16 iterations indepen- 
dent of the resolution N . 

How much computational work is required per iteration? 

The righthand side of (6.4) can be evaluated in 0 (N tog N) 
operations by transform methods because u^ n ^ is known. Also 
ithe solution of (6.4) for u^ n+1 ^ can be accomplished 
^efficiently because of assumption (ii) above. 

Chebyshev iteration 

The rate of convergence of the scheme (6.4) can be 
accelerated using Chebyshev acceleration [8] . The scheme is 


(n+1 ) 


ap N 


L [oj u< n) 
ap n N 


+ ( 1 — U) 


) u^ n 1 ^]-aaj (L u^. n ^ 
N n sp N 


-V 


( 6 . 10 ) 


where 


0) 

n 


23T n (6) 

T n + l (e) 


( 6 . 11 ) 


and 8 = min(|l-am| 1 , |l-aM| ^ ) . Here T n (8) is the 
Chebyshev polynomial of degree n and m,M are given in 
(6.3) . It is not hard to show that the error in the solution 
of (6.1) decreases by at least the factor 8 + /p-1 after 
each iteration of (6.10) . 

If M = 2.5, m = 1, then choosing a= 4/7 gives 
8 = 7/3, so the error decreases by at least ( 7 + 2/TO) /3~4 . 44 
at each iteration. Therefore the error in the solution of 
(6.1) is decreased by a factor 10^ after about 9 iterations 
of (6.10) , nearly a factor two faster than the Richardson 
method (6.4). The penalty of using the Chebyshev acceleration 
method is that two levels of iterates, u^ n ) and 
must be stored. 
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Conjugate gradient iterations 


For many problems, it is possible to accelerate the 

convergence of the iteration method for solving ( 6 . 1 ) still 

further by using the conjugate gradient method [ 9 ] . The 

applicability of this method to the solution of spectral 

equations will be studied in depth elsewhere. The result is 

that for general Dirichlet problems for elliptic partial 

4 

differential equations, the error is decreased by 10 after 

7 

only 3 conjugate gradient iterations and by 10 after 
only seven iterations. 

Methods to construct a suitable approximation operator 

L ^ will now be explained. The idea is akin to an idea of 

ap 

D'yakonov [10], but differs in one essential respect. 
D'yakonov considers the solution of the elliptic partial 
differential equation 

Du = V* K (x , y ) V u = f (6.12) 

in some region R with, say, Dirichlet boundary conditions 
on <5 r. If 


< K . < K (x , y) < K 

mm — J — max 


then D'yakonov proposes the iteration scheme 


2 (n+ 1 ) 


u 


2 (n) 
= V u 


a (Du 


(n) 


- f) 


(6.13). 


where a- 2/ (K . +K ) . The error ||u^ n ^-u II decreases 
mm max 11 11 

by at least the factor (K -K . ) / (K +K . ) at each 

max min ' max min 

iteration of (6.13). D'yakonov 1 s method involves approx- 
imating the differential operator by a new, simpler, operator 
whose spectrum is bounded in terms of the original operator. 

L is constructed from L by changing 
ap sp 2 2 2 

the discretization operator either in addition to or in place 

of approximating the differential operator. Thus, L 
is constructed by a suitable low-order finite differenc e app- 
roximation to L. 
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A simple example is given by the second-order differential 
equation 

Lu = f(x)u" (x) + gCx)u'(x) + h(x) = v (x) (0<x<2tt) 

(6.14) 

with periodic boundary conditions u(x+2tt) = u(x) and 
f (x) >0. A spectral approximation is approximately sought as 
the finite Fourier series 


u ( x ) = Y a, e ikx . (6.15) 

| k I <K K 

If the Fourier coefficients of f (x) , g(x), h(x), v (x) are 
denoted f^> 9 ]^ ' ' V k ' respectively, then the spectral 
(Galerkin) equations for are 


L u = 
sp 


,1 

p|<K 
k-p < K 


2 _ 
■P f 


- n +i P9lc-, 


+h. 


] a = v 


k-p ’ "^k-p ’ **k-p J “p k’ 


(6.16) 


Clearly, these equations have, in general, a full matrix 

2 

representation that requires 0(K ) storage locations and 
O(K^) operations to invert. 

A suitable approximate operator L construc ‘* ::ec ^ 

using the collocation points x^ = 2-rrj/N ( j = 0 , 1 , . . . , N-l ) 
where N = 2K. In the physical space representation, 
the finite difference approximation 


L u 
ap 


f (Xj) 


U • . -| -2u . H-U . n 

_I±1 2 lzl 

(Ax) 2 


+g(x j ) 


u.,..-u. 

. j +1 ml 

2 Ax 


+h (Xj ) 


(6.17) 


and Ax = 2 tt/N is used. 


where u . = u (x . ) 

3 2 

sparse and efficiently invertible. To verify 


Obviously, L is 
J ap 

(6.3) we use 


the following elementary argument (that may be made more 
rigorous but no more correct when more involved WKB-like 
arguments are used.) If X is an eigenvalue of L then 

there exists a function u(x) such that 


L u = X L u (6.18) 

sp ap 
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If u(x) is a smooth function of x (in the limit N , 
then both L u and L u should be good approximations 

to Lu(x) so (6.18) implies x~l- On the other hand, if 
u (x) is a highly oscillatory function of x (in the limit 
N 00 ) then 


u" >> u ' >> u (N -»- °°) . 


(6.19) 


Therefore , 


L u ~ f (x . ) 
ap j 


u i+ 1 - 2u i~ fu i _ 1 

(Ax) 2 


( 6 . 20 ) 


and, if transform (pseudospectral) methods are used to eval- 
uate L u, 
sp 


_ . . , , 2 » ikx ■ 

L u ~ f (x . ) ) (-k ) a, e j 

sp 3 |)c|<K k 


( 6 . 21 ) 


so (6.18) gives 


f (x . ) 




k <K 


, ,2. rkx • . _ , . 

(-k ) a, e Xf(x. ) 

J 


u, ,.-2u.+u. .. 
3 +1 3 1^1 

(Ax) 2 


( 6 . 22 ) 


The eigenfunctions of (6.22) are 

Uj = e^^j ( j q | <K ) (6.23) 

and the associated eigenvalue is 

X = (6.24) 

4 sin 2 q Ax 

Since |q|<K with K = ^ N = tt/Ax, it follows that 
2 

1< X < \ (6.25) 

2 

Thus, (6.3) holds with m = 1 and M = rr /4 ~ 2.5. 

There are several extensions of the above method for 
constructing L that are important in practice. 


First, in 



1 


the case of Chebyshev-spectral methods, it is appropriate to 
construct L using finite-difference approximations based 
on the collocation points x^ = cos irj/N. In this case, the 
operator bounds (6.3) continue to hold with M = 2.5, m = 
for a wide variety of operators L. Second, higher order 
equations are best treated by writing them as a system of 


lower order equations. Thus, direct construction of L 


ap 


L = V 


gives 


for 


X < | |L 


-1 

ap 


2 2 


sp 


< 6 


( 6 . 26 ) 


2 

However, introducing v = V u, defining the second-order 
operator K by 


K( U ) 

v 


V 

V 


2 

u-v 


(6.27) 


and direct construction of 
operator gives 



1 i I |K 


-1 

ap 



<2.5. 


as a finite-difference 


(6.28) 


Third, odd-order operators, initial-value problems, and prob- 
lems of mixed type are best treated by constructing L on 
a grid that is roughly 50% finer than that used in construction of 
L by collocation. In this case the spectral bounds (6.3) 

with M< 2.5 continue to hold for most problems. For example, 

3 

the operator with periodic boundary conditions has spectrum 

ik while its centered finite-difference approximation has 
spectrum i sin(kAx)/Ax so 


I |l ^ L 11= 0(kAx/sin kAx) 

' 1 ap sp 1 

which is unbounded for jkAx|<7T, but bounded by 4 it/ 3/3~2.4 
if |kAxj < 2ir/3. 

7. EXAMPLES 

Consider the solution of the heat equation 

= y 2 u - e x - e y (x,y)eD (7.1) 
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with Dirichlet boundary conditions 


u(x,y) = e x + e y (x,y) s3D, (7.2) 

where D is the annular region plotted in Fig. 1 whose inner 
and outer boundaries are, respectively, (in polar coordinates) 

r = f^(0) = 0.3 + 0.1 sin 0 + 0.15 sin 50 (7.3) 


and 


r = f 2 ( 0 ) =1 + 0. 2 co S0+0. 15s in 4 0. (7.4) 

This problem is solved spectrally using the stretching trans- 
formation (5.3) to transform D into the rectangular domain 
(5.2) and then using Fourier series in 0 and Chebyshev 
series in z to represent u(x,y) . 

As t -*■ “, the solution to (7.1) - (7.4) approaches the 
steady state solution 


X V 

u(x,y) = e + e J (x,y)e D 


(7.5) 


In Table 1, the maximum pointwise errors in the spectral 
solution of this problem as t -> 03 are given. Evidently the 
rapid convergence of the spectral solution to the exact steady 
state is preserved in this complex geometry problem. It is 
not very surprising that, for a given total number of retained 
modes, the best accuracy is achieved with about 4 times more 
'angular 1 (0) than 'radial' (z) modes; after all, the 
'annulus' (7.3) - (7.4) is on average about tt times longer 
in circumference than in radius. 

Another example of the technique suggested in Sec. 6 is 
the solution of the boundary-layer equations 


~ + f + 3(0 [i - <!4 )2] 
an 3 3n Sn 


ril d ^l - ill ill 
^ l 3n 3C3n g 2 3£ J 

(7.6) 


where f(C,n) is the dimensionless stream function and is sub- 
ject to the boundary conditions 



f (5/0) 


(7.7) 


- H < ? ' 0) “ ° 

~ U,n> + 1 (n -*■ °°) (7.8) 

3n 

Here ( ^ , n ) are the Levy-Lees coordinates; E, increases in the 
free-stream direction and n increases away from the wall. The 
non-self -similar solution of these equations for Howarth’s 
flow in which the pressure-gradient parameter (3(£) is given 
by 


6(5) = ^4 (7.9) 

is a standard test problem [11]. Eqs. (7.6) - (7.9) are 

solved using Chebyshev series in the transformed variable 

s = ^ - i <-i < s < i) , 

where R is chosen so r\ = R is in the free stream and a Crank- 
Nicolson scheme is used in £ . The finite-difference errors 
in 5 are reduced by Richardson extrapolation [13] . The 
results for the viscous wall stress at two downstream locat- 
ions near the separation point £ a 0.901 are given in Table 
2 together with some corresponding finite difference results 
from Ref. 13. These results illustrate the high accuracy 
of the spectral schemes with modest resolution. 

8 . OTHER DEVELOPMENTS 

In this paper, techniqueshave been explained to extend 
spectral methods efficiently to problems in complex geometries. 
Many applications of these techniques are underway and results 
will be presented elsewhere. 

There have been several other recent developments of 
spectral methods that relate to these applications but have 
independent interest. First, Lagrangian spectral methods 
have been developed for the solution of high speed flows. In 
this case, Lagrangian marker particles are used to provide the 
coordinate transformations necessary for the methods introduced 
in this paper. Second, fast conformal mapping techniques have 
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been developed. These techniques are particularly attractive 
for the solution of free-surface flow problems by spectral 
methods. 
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Table 1. Errors in Steady State Solution of (7.1) -(7. 4) 


Number of 

Number 

Maximum 

angular 

of ’radial ' 

relative 

modes (0) 

Chebyshev modes (z) 

error 

16 

4 

1.4 x 10" 2 

32 

8 

2.8 x lcT 5 

64 

16 

2.5 x 10“ 10 
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Table 2. viscous Wall Stresses for Howarth's Flow 


Method Truncation 

Numbdr of grid 

Number of grid 

Total 

A 

B 

R 

0<_s<R 

points (modes) 

points in 

number 



Method 


in n 

i 

of points 



Spectral 

8 

23 

6 

138 

0.21239 

0.04454 

Spectral 

8 

23 

11 

253 

0.20913 

0.03471 

Spectral with 
Richardson 







extrapolation 
in £ 

8 

23 

{6,11} 

253 

0.20718 

0.03064 

Finite 

difference 

13] 

6 

121 

51 

6171 

0.20723 

0.03083 

Finite 
dif f erence 







with Richardson 






extrapolation 
in £ and/or n 
[13] 

6 

19 

16 

489 

0.20710 

0.03171 


6 

61 

16 

1708 

0.20741 

0.03121 


6 

121 

51 

9282 

0.20714 

0.03053 


A = <£= 0.7, n = 0) ~ 0.20724 B = (£= 0.894,n = 0) - 0.03072 

3n 3n 
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